#include <assert.h>
#include <stdio.h>
#include <stdlib.h>

#include "../src/em.c"

typedef struct _EmTestDatum EmTestDatum;
struct _EmTestDatum 
{ 
  struct
    {
      float     kappa;
      float     logalpha;
      float     tau[(5-1) * 4];
      float     loglambda[4];
      float     priorz[5];
      int       clamp;
      Rating    ratings[8];
      float     logalphastddev;
    }                   inputs;

  struct
    {
      float    newq;
      float    logalpha;
      float    pz[5];
    }                   desired_output;
};

static void
test_em (void)
{
  EmTestDatum data[] = 
    {{{0.923811,  0.498408, {0.0112682, 0.634888, 0.460961, 0.110676, 0.467752,  0.947949, 0.821227, 0.145074, 0.257316, 0.862082, 0.161289,  0.281663, 0.517215, 0.0815973, 0.96935, 0.634766}, {0.698277,  0.283516, 0.226767, 0.743543}, {0.893733, 0.655741, 0.302956,  0.245136,  0.882464}, -1, {{3, 0}, {3, 2}, {1, 0}, {0, 2}, {3, 3}, {1,  1}, {3, 2}, {2, 1}},  1.842}, {-12.669294912742228159, -1.0697869669947522034,  {0.0163861, 0.432428, 0.513334, 0.0374504, 0.000401675}}}, {{0.13446,  0.414712, {0.072903, 0.0207682, 0.989386, 0.157396, 0.210821,  0.859479, 0.707723, 0.640181, 0.129223, 0.890129, 0.0729568,  0.941904, 0.845707, 0.663362, 0.329413, 0.0481713}, {0.189966,  0.360405, 0.0842774, 0.165707}, {0.169114, 0.51841, 0.949818,  0.750995,  0.0962107}, -1, {{0, 0}, {0, 2}, {0, 2}, {3, 4}, {0, 3}, {2,  2}, {0, 0}, {0, 3}},  1.96043}, {-13.944306521052668895, -1.7458488915506438603,  {0.0177388, 0.153546, 0.464616, 0.341204, 0.022895}}}, {{0.5936,  0.88539, {0.638163, 0.252709, 0.953419, 0.756167, 0.748034,  0.179752, 0.0115151, 0.91046, 0.084672, 0.850339, 0.963344,  0.720494, 0.724266, 0.766061, 0.797637, 0.55138}, {0.205856,  0.816244, 0.0466414, 0.45517}, {0.708215, 0.855812, 0.453042,  0.569779, 0.0700519},  1, {{2, 3}, {2, 2}, {3, 0}, {2, 1}, {1, 4}, {2, 1}, {3, 0}, {1,  2}}, 1.49962}, {-15.662050283412319642, -1.8205865261973466137,  {8.39107e-7, 0.999992, 3.64457e-6, 2.38157e-6,  7.06254e-7}}}, {{0.813612,  0.322018, {0.423352, 0.488108, 0.903152, 0.237346, 0.573013,  0.524764, 0.182658, 0.51308, 0.806952, 0.727127, 0.631278,  0.307224, 0.990708, 0.680486, 0.176108, 0.599009}, {0.134896,  0.227444, 0.606329, 0.528957}, {0.531793, 0.727821, 0.792716,  0.206939,  0.108441}, -1, {{3, 0}, {3, 0}, {1, 4}, {1, 0}, {3, 4}, {0,  2}, {2, 3}, {2, 3}},  1.88956}, {-15.195218433699467574, -2.4304252378035908298,  {0.114955, 0.286176, 0.432463, 0.117892, 0.0485136}}}, {{0.969592,  0.535428, {0.714949, 0.706906, 0.456512, 0.728476, 0.987822,  0.0756282, 0.149289, 0.737768, 0.307336, 0.89952, 0.55028,  0.602871, 0.0798915, 0.293191, 0.0213225, 0.0710787}, {0.35207,  0.500475, 0.814384, 0.962638}, {0.112357, 0.610911, 0.844792,  0.42721, 0.397407},  1, {{2, 1}, {0, 4}, {0, 1}, {0, 3}, {0, 2}, {3, 3}, {3, 2}, {3,  0}}, 1.38828}, {-16.681424503748581042, -1.8846996896412207069,  {4.47765e-7, 0.999991, 4.93508e-6, 2.98139e-6,  5.74641e-7}}}, {{0.698734,  0.409585, {0.828377, 0.238991, 0.960966, 0.10225, 0.928857,  0.688711, 0.358095, 0.0223582, 0.635665, 0.667388, 0.287016,  0.670288, 0.13519, 0.853005, 0.324379, 0.557931}, {0.524279,  0.0082131, 0.897169, 0.160524}, {0.620275, 0.619934, 0.198435,  0.750938,  0.791898}, -1, {{0, 2}, {1, 0}, {0, 3}, {2, 4}, {1, 1}, {0,  1}, {2, 4}, {0, 0}},  1.23747}, {-16.043887465279793278, -1.8031552958438401694,  {0.00316146, 0.0434938, 0.0759775, 0.530352, 0.347015}}}};

  for (unsigned int i = 0; i < sizeof (data) / sizeof (data[0]); ++i)
    {
      ParameterizedNormalDistribution prior_logalpha = 
        parameterized_normal_distribution (data[i].inputs.logalphastddev);

      float logalpha = data[i].inputs.logalpha;
      float pz[5];
      float logpriorz[5];

      for (unsigned int j = 0; j < 5; ++j)
        {
          pz[j] = drand48 ();
          logpriorz[j] = log (data[i].inputs.priorz[j]);
        }

      float q = em (data[i].inputs.kappa,
                    &logalpha,
                    5,
                    data[i].inputs.tau,
                    data[i].inputs.loglambda,
                    pz,
                    logpriorz,
                    data[i].inputs.clamp,
                    data[i].inputs.ratings,
                    8,
                    &prior_logalpha.base,
                    1000,
                    1e-4);

      assert (fabs (q - data[i].desired_output.newq) <=
              1e-3 * (1 + fabs (q) + fabs (data[i].desired_output.newq)) ||
              ( fprintf (stderr, "%u: %g ?= %g\n", i, q, data[i].desired_output.newq),
                0 ));
      
      assert (fabs (logalpha - data[i].desired_output.logalpha) <=
              1e-3 * (1 + fabs (logalpha) + fabs (data[i].desired_output.logalpha)) ||
              ( fprintf (stderr, "(oops logalpha ... but this is ok) %g ?= %g\n", logalpha, data[i].desired_output.logalpha),
                1 ));

      for (unsigned int j = 0; j < 5; ++j)
        {
          assert (fabs (pz[j] - data[i].desired_output.pz[j]) <=
                  1e-3 * (1 + fabs (pz[j]) + fabs (data[i].desired_output.pz[j])) ||
                  ( fprintf (stderr, "(oops pz[%u] ... but this is ok) %g ?= %g\n", j, pz[j], data[i].desired_output.pz[j]),
                    1 ));
        }
    }
}

int 
main (void)
{
  srand48 (69);

  test_em ();

  return 0;
}
